rm(list = ls())

setwd('/path/to/replication/')

library(data.table)
library(ggplot2)

# Figure B4: Predicted vs. observed probability of death
excess_deaths_all <- get(load(file='./data/excess_deaths_cdc_methodology.RData'))
ps <- glm(excess_deaths_feb_d ~ excess_deaths_jan, family=binomial, data=excess_deaths_all)
excess_deaths_all$prob_excess_deaths <- predict(ps, type="response")

excess_deaths_all[prob_excess_deaths<.1,bin:=.1]
excess_deaths_all[prob_excess_deaths>=.1 & prob_excess_deaths<.2,bin:=.2]
excess_deaths_all[prob_excess_deaths>=.2 & prob_excess_deaths<.3,bin:=.3]
excess_deaths_all[prob_excess_deaths>=.3 & prob_excess_deaths<.4,bin:=.4]
excess_deaths_all[prob_excess_deaths>=.4 & prob_excess_deaths<.5,bin:=.5]
excess_deaths_all[prob_excess_deaths>=.5 & prob_excess_deaths<.6,bin:=.6]
excess_deaths_all[prob_excess_deaths>=.6 & prob_excess_deaths<.7,bin:=.7]
excess_deaths_all[prob_excess_deaths>=.7 & prob_excess_deaths<.8,bin:=.8]
excess_deaths_all[prob_excess_deaths>=.8 & prob_excess_deaths<.9,bin:=.9]
excess_deaths_all[prob_excess_deaths>=.9,bin:=1]

excess_deaths_all[,c('p','n'):=.(mean(excess_deaths_feb_d),.N), by=bin]
calibration_plot <- unique(excess_deaths_all[, .(bin, p, n)])
calibration_plot[,sd:=sqrt((p*(1-p))/n)]
calibration_plot[,c('ci.low', 'ci.high'):=.(p-1.96*sd, p+1.96*sd)]

ggplot(calibration_plot, aes(x=bin, y=p)) + 
  geom_point(colour='#1b2631') +
  geom_line(colour='#154360') +
  geom_ribbon(aes(ymin = ci.low, ymax = ci.high), alpha=0.25) +
  scale_x_continuous(breaks=seq(-.25,1,.25),limits = c(-.25,1)) +
  scale_y_continuous(breaks=seq(-.25,1,.25),limits=c(-.25,1)) +
  geom_abline(intercept = 0, slope = 1, linetype='dashed',colour='#666666') +
  ylab('Observed probability of health threat') +
  xlab('Predicted probability of health threat') +
  theme(
    legend.position="none",
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black"
    ),axis.title.y = element_text(size = rel(1.2)), axis.title.x = element_text(size = rel(1.2)))
ggsave('./reliability_curve_health.pdf', width=6, height=4.5)


# Figure B5: Average number of excess deaths by month
death_series <- fread('./data/municipality_monthly_deaths_2007_2020.csv')

death_series[order(panelvar, month), deaths_lag12m:=shift(deaths, n=12, type="lag"), by=.(panelvar)]
death_series[, deaths_excess:=deaths-deaths_lag12m]
death_series[order(panelvar, month), deaths_excess_lag1m:=shift(deaths_excess, n=1, type="lag"), by=.(panelvar)]

ggplot(
  death_series[!is.na(deaths_excess) & month<'2020-03-01', mean(deaths_excess), by=month]
) + aes(x=month, y=V1) + geom_point() + geom_line() + ylab('average municipality excess deaths') +
  xlab(NULL) + geom_vline(xintercept = as.Date('2020-01-01'), linetype='longdash', color='grey') +
  theme(
    legend.position="none",
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black"
    ),axis.title.y = element_text(size = rel(1.1)), axis.title.x = element_text(size = rel(1.2)))
ggsave('./excess_deaths_by_month.pdf', width=6, height=4.5)


# Figure B6: Predicted vs. observed probability of unemployment
load('./data/panel_month_dummies.RData')
unemployment <- fread('./data/unemployment_region.csv')
treat_cities <- unique(panel[,.(panelvar, region, losing_income_d)])
propensity <- merge(treat_cities,unemployment,by='region',all.x = T)
ps <- glm(losing_income_d ~ unemployment, family=binomial, data=propensity)
propensity$prob_losing_income <- predict(ps, type="response")

propensity[prob_losing_income<.1,bin:=.1]
propensity[prob_losing_income>=.1 & prob_losing_income<.2,bin:=.2]
propensity[prob_losing_income>=.2 & prob_losing_income<.3,bin:=.3]
propensity[prob_losing_income>=.3 & prob_losing_income<.4,bin:=.4]
propensity[prob_losing_income>=.4 & prob_losing_income<.5,bin:=.5]
propensity[prob_losing_income>=.5 & prob_losing_income<.6,bin:=.6]
propensity[prob_losing_income>=.6 & prob_losing_income<.7,bin:=.7]
propensity[prob_losing_income>=.7 & prob_losing_income<.8,bin:=.8]
propensity[prob_losing_income>=.8 & prob_losing_income<.9,bin:=.9]
propensity[prob_losing_income>=.9,bin:=1]

propensity[,c('p','n'):=.(mean(losing_income_d),.N), by=bin]
calibration_plot <- unique(propensity[, .(bin, p, n)])
calibration_plot[,sd:=sqrt((p*(1-p))/n)]
calibration_plot[,c('ci.low', 'ci.high'):=.(p-1.96*sd, p+1.96*sd)]

ggplot(calibration_plot, aes(x=bin, y=p)) + 
  geom_point(colour='#1b2631') +
  geom_line(colour='#154360') +
  geom_ribbon(aes(ymin = ci.low, ymax = ci.high), alpha=0.25) +
  scale_x_continuous(breaks=seq(0,1,.1),limits = c(0,1)) +
  scale_y_continuous(breaks=seq(0,1,.1),limits=c(0,1)) +
  geom_abline(intercept = 0, slope = 1, linetype='dashed',colour='#666666') +
  ylab('Observed probability of unemployment threat') +
  xlab('Predicted probability of unemployment threat') +
  theme(
    legend.position="none",
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black"
    ),axis.title.y = element_text(size = rel(1.2)), axis.title.x = element_text(size = rel(1.2)))
ggsave('./reliability_curve.pdf', width=6, height=4.5)

# Figure B7: Correlation between local exposure to death and unemployment
load('./data/panel_month_dummies.RData')
panelu <- panel[month=='2019-12-01', .(panelvar, excess_deaths_jan_d, losing_income_d)]
panelu <- panelu[!is.na(excess_deaths_jan_d)]

df = structure(list(deaths = structure(c(1L, 1L,2L, 2L), .Label = c("high deaths", 
                                                                    "low deaths"), class = "factor"), unemployment = structure(c(1L, 2L, 1L, 2L), 
                                                                                                                               .Label = c("high unemployment",  "low unemployment"), class = "factor"),
                    n = unlist(c(panelu[,.N, by=.(excess_deaths_jan_d, losing_income_d)][1,3],
                                 panelu[,.N, by=.(excess_deaths_jan_d, losing_income_d)][4,3],
                                 panelu[,.N, by=.(excess_deaths_jan_d, losing_income_d)][3,3],
                                 panelu[,.N, by=.(excess_deaths_jan_d, losing_income_d)][2,3]))), 
               row.names = c(NA, -4L), class = "data.frame")

xtab <- xtabs(n~unemployment+deaths, data=df)
pdf('./cor_death_unemployment.pdf',width=5, height=5.3)
mosaicplot(xtab, col=TRUE, main=NULL, cex.axis = 1)
dev.off()

# Figure B10: Correlation between vote for the extreme right in national elections and far-right mayors
### correlation between far-right mayors and extreme right vote share
load('./data/panel_month_dummies.RData')
panelu <- panel[month=='2019-12-01', .(panelvar, extr_right_mayor, extr_right_d)]
panelu <- panelu[!is.na(extr_right_mayor) & !is.na(extr_right_d)]

df = structure(list(mayor = structure(c(1L, 1L,2L, 2L),
                                      .Label = c(paste("with far-right", "mayor", sep='\n'), 
                                                 paste("without far-right", "mayor", sep='\n')), class = "factor"), voteshare = structure(c(1L, 2L, 1L, 2L), 
                                                                                                                                          .Label = c(paste("high extreme", "right vote share", sep = '\n'),
                                                                                                                                                     paste("low extreme", "right vote share", sep='\n')), class = "factor"),
                    n = unlist(c(panelu[,.N, by=.(extr_right_mayor, extr_right_d)][4,3],
                                 panelu[,.N, by=.(extr_right_mayor, extr_right_d)][3,3],
                                 panelu[,.N, by=.(extr_right_mayor, extr_right_d)][2,3],
                                 panelu[,.N, by=.(extr_right_mayor, extr_right_d)][1,3]))), 
               row.names = c(NA, -4L), class = "data.frame")

xtab <- xtabs(n~voteshare+mayor, data=df)
pdf('./cor_mayor_voteshare.pdf',width=5, height=5.3)
mosaicplot(xtab, col=TRUE, main=NULL, xlab = 'vote share', cex.axis = 1, las = 2)
dev.off()
